Self-organized criticality in a rice-pile model 
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We present a new model for relaxations in piles of granular material. The relaxations are 
determined by a stochastic rule which models the effect of friction between the grains. We find power- 
law distributions for avalanche sizes and lifetimes characterized by the exponents r = 1.53 ± 0.05 
and y = 1.84 ± 0.05, respectively. For the discharge events, we find a characteristic size that scales 
with the system size as , with [i — 1.20 ± 0.05. We also find that the frequency of the discharge 
events decrease with the system size as with /i' — 1.20 ± 0.05. 

PACS numbers: 05.40. +j, 64.60. Ht, 05.70.Jk, 05.70.Ln 



Since its introduction by Bak, Tang, and Wiesenfeld 
0, the concept of self-organized criticality (SOC) and 
models which display SOC behavior have been the fo- 
cus of much interest However, comparison with 
real systems has proved to be a though test for the the- 
ory and models ||lo|-[l^. Furthermore, in one dimension 
the models tend to display either trivial behavior or be- 
havior that cannot be classified as critical. Against this 
background, recent experiments on rice piles ]l^ have 
shown that under some conditions a real rice pile can 
self-organize into a critical state: For grains with a large 
aspect ratio the system self-organizes into a critical state. 
Frette et al. explained this result with the increased fric- 
tion and packing possibilities that were able to cancel 
inertia effects. Furthermore, they observed that large lo- 
cal slopes developed in the pile. 

Here, we propose a new model for a pile of granular 
material where we introduce randomness in the relax- 
ation rule instead of in the deposition rule. We study the 
model in one dimension and find power law distributions 
for avalanche sizes s and lifetimes T. We also study the 
distribution of sizes for discharge events (i.e., particles 
falling off the pile), and find it to be bounded. The re- 
sults show that our model belongs to a new universality 
class for systems displaying SOC. 

First, we define the one-dimensional model. The sys- 
tem consists of a plate of length L, with a wall at i = 
and an open boundary at i = L + I. The profile of the 
pile evolves through two mechanisms: deposition and re- 
laxation. Deposition is always done at i = 1, and one 
grain at a time. The rate of deposition is slow enough 
that any avalanche, initiated by a deposited grain, will 
have ended before a new grain is deposited. 

During relaxation we look at all active columns of the 
rice pile: A column i of the pile is considered active if, in 
the anterior time step, it (i) received a grain from column 
« — 1, (ii) toppled a grain to column i -I- 1, or (iii) column 
i + l toppled one grain to its right neighbor. If a column i 
is active and the local slope, i.e., 6h{i) = h{i) — h{i + 1), 
is strictly larger than a threshold value 5*1, then with 
probability p a grain will move from i to i + 1. However, 
if Sh{i) > 5*2, a grain is moved from i to i + 1 with prob- 
ability one. Grains toppled from column i — L leave the 
system. When no active columns remain on the pile, the 
avalanche is said to be over. 
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FIG. 1. a) In the inset, we show the average slope as 
a function of time (measured as the number of deposited 
grains). As all the data presented in this paper, the aver- 
age slope was obtained for p = 0.6, = 1 and 5*2 = 4. The 
system sizes shown are L = 40, 80, 160. We can see that after 
an initial transient regime, whose duration depends strongly 
on the system size, a steady state is reached. The figure shows 
the fluctuations in the average slope of the pile, and the size 
of the discharge events, as a function of time for as system of 
size 80. It is visually apparent that there is a close connection 
between sharp changes in the average slope and the discharge 
events at the boundary. It is interesting to note that although 
our model does not include either a repose or a maximum an- 
gle for the pile, the dynamics suggest the existence of such 
angles because of the large discharge events, b) Plot of the 
avalanches sizes as a function of time for the same system and 
interval as in (a). Again, the connection between the largest 
avalanches and the discharge events is observed. 
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The physical interpretation of our rules is the follow- 
ing: Suppose that a column, or portion, of the pile is in a 
metastable configuration. If a new grain is deposited or 
toppled on top of it, or the local slope changes, then that 
metastable configuration can become unstable. To model 
such an effect, we introduce the parameter p, which repre- 
sents the fact that there is a finite probability that a new 
stable configuration is reached. Physically the parameter 
p thus describes the friction between the rice grains and 
the possibility that a metastable packing configuration 
will be attained. The friction effect is the main new in- 
gredient in our model compared to other models and it 
comes directly from the observation that there exists a 
large range of slopes in the rice pile instead of a single 
critical value The friction p can be a complicated 

function of local slopes and packings of the particles but 
we find that the results are insensitive to the specific 
form and value of p. The parameter 5*2 models the ef- 
fect of gravity on the packing arrangements. We assume 
that above the maximum value S2 of the local slope, it is 
no longer possible for a local stable configuration to be 
achieved, thus a grain must be toppled. In the limiting 
cases p = 0, 1, or S'2 = ^i, we recover the model in Ref. 
(which has trivial behavior for one dimension). 

The simulation of the model shows two distinct 
regimes, a transient period followed by a steady (criti- 
cal) state; cf. Fig. [l]. Here, we focus on the properties 
of the model in the critical state. As can be seen in 
Fig. 1^, the model leads to the establishment of a state 
with wildly varying avalanches sizes and a complicated 
structure in time. The size of an avalanche can be de- 
fined in a number of ways: the number of topplings s, the 
lifetime of the avalanche T, or the size of the discharge 
events m. We start by investigating the distribution of 
s. Figure ^(a) shows the probability density of avalanche 
sizes for different system sizes. The distribution follows 
the scaling form 



rid ~ L~ 



(3) 



P{s,L) 



fs{s/Ln, 



(1) 



where fs is a scaling function rapidly decaying for large 
arguments. The best collapse is obtained with the ex- 
ponents T = 1.53 ± 0.05 and i/ = 2.20 ± 0.05, cf. 
(b). Even though r is close to the mean-field value 
15|-|l7[|, our model describes a new universality class. 
This can be shown by mapping the avalanche dynamics to 
the motion of an interface through a disordered medium 
[p^ . By using that the average number of topplings is 
(s) = L in the critical state, it follows from (|l|) that 




r = 2-i. 



(2) 



in agreement with our numerical results. 

An interesting characteristic of the distribution is the 
presence of a peak, deviating from the power-law behav- 
ior, for a size close to the cutoff of the distribution. A 
close look at Fig. |l| shows that the biggest avalanches 
coincide with large changes in the average slope of the 
pile and with discharge events. Furthermore, as shown 
in Fig. 1^, the number Ud of avalanches reaching the open 
boundary, for a given number of deposited grains, scales 
with the system size as 



with /i' ~ 1.20 ±0.05. This suggests that the peak is due 
to finite-size effects which lead the system into a super- 
critical state, followed by a massive avalanche and a large 
change in the average slope. We check this hypothesis by 
considering only the avalanches for which no discharge 
occurred at the boundary. As can be seen in Fig. |^(a) , the 
peak for very large values of s is then no longer present 
and the cutoff has moved to a smaller value, confirming 
our hypothesis. We find that the data is described by 
Eq. (0) with the same values of the exponents as when 
the peak is present. 
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FIG. 2. a) Log-log plot of the probability density of 
avalanche sizes s for several system sizes. It is visually ap- 
parent that for s 1 a power-law dependence is observed. 
For values of s close to the cutoff, imposed by the finiteness of 
the system, we observe a peak deviating from the power-law 
behavior. To show that the peak is due to the system reach- 
ing a supercritical state (which is followed by a discharge) 
we also plot the distribution of avalanches when no discharge 
event occurred: That curve, obtained for L = 640, is shifted 
vertically by a factor of 8, to make it more visible, b) Data 
collapse of the curves shown in (a) according to Eq. (nl) with 
the exponents r ~ 1.53 and v ~ 2.20. 

Next, we study the distribution of lifetimes T for the 
avalanches. As shown in Fig. ^(a), the data is described 
by the scaling form 



P{T,L)^T-y /t(T/L'^), 



(4) 



which is confirmed by the good data collapse obtained 
with the exponents y = 1.84 ± 0.05 and a = 1.40 ± 0.05. 
From conservation of probability follows that a{y — 1) = 



2 



i/(t — 1) ||T^, in nice agreement with our results. Finally, 
we study the distribution of sizes m for the discharge 
events (see Fig. |](a)). The scaling ansatz 



(5) 



where the scaling function /,„ decays exponentially, de- 
scribes the data. Since the distribution P(m, L) does not 
diverge for m — > and its integral must equal one, it fol- 
lows that K = fi. This result is confirmed by the data 
collapse shown in Fig. ||(b), obtained for the exponents 
K = 1.2 ±0.1 and n = 1.2 ±0.1. 




FIG. 3. The dependence of the number of discharge events 
Ud as a function of L. A power law behavior, cf. Eq. (^), with 
ji' ~ 1.20 is obtained. 

It is possible to obtain additional scaling relations for 
the exponents besides those mentioned above. In the 
steady state, the input of matter must balance the out- 
put through the open boundary. Thus, we obtain that 
the frequency of discharge events must balance their char- 
acteristic size, and fi = fj,' . The characteristic size of the 
discharge events depends on the system size as L^. So, 
we can conclude that whenever the system reaches a su- 
percritical state, the number of grains discharged is of 
order L^^. Since the average number of topplings for a 
given grain before being discharged is of order L, it fol- 
lows that the cutoff size for the avalanches must scale as 
L X ^ L'^, thus 1/ = 1 + fi, in accordance with our 
results. 
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FIG. 4. a) Log-log plot of the probability density of 
avalanche lifetimes T for several system sizes. It is visually 
apparent that for t S> 1 a power-law dependence is observed. 
As for the avalanche sizes, a peak is present for lifetimes close 
to the cutoff, b) Data collapse of the curves shown in (a) ac- 
cording to Eq. (^) with the exponents y ~ 1.84 and a ~ 1.40. 

Just before sumbission of the present work, we became 
aware of a model by Christensen et al. |l9| which for some 
range of parameters seems to belong to the same univer- 
sality class as the model discussed here. The model in 
Ref. introduces stochasticity in the toppling of par- 
ticles via the selection of a new random critical slope 
for columns where a toppling occurred. For local slopes 
above the critical slope, a grain is always toppled. In 
Ref. |l9) the predictions of the model are compared with 
experimental results for the diffusion of tracer particles. 
The numerical value of the exponent describing the diffu- 
sion of the tracers is in rough agreement with the numer- 
ical predictions of the model. The exponent a describing 
the scaling of the potential energy dissipated during an 
avalanche is a ~ 1.53 for our model pO| , in disagree- 
ment with the experimental result a ~ 2. In fact, the 
large disagreement between the experimental value of a 
and the numerical prediction of the model suggest that 
our model and the model of Ref. |^ do not belong to 
the same universality class as the rice-pile experiment. 
In a stochastic sandpile model is studied in which 
only the front of the avalanches propagate (i.e., no back- 
ward avalanches are allowed) . Such a rule leads to higher 
values for r but apparently at the cost of destroying uni- 
versality. 

Frette et al. also found that for "round" rice grains the 
system did not evolve into a critical state, and that the 
distribution of avalanche sizes was bounded. The reason 
for this result can be understood if the results for the role 
of inertia on the dynamics of sandpiles are remembered 
1^,^. It was shown in Ref. Q that the assumption 
of zero inertia is essential for the establishment of the 
critical state. For real experiments, where inertia cannot 
be avoided, that assumption can only be valid for sys - 
tem with sizes smaller than a threshold value Lc p2| , p3| . 
Thus, for the round rice grains all system sizes studied 
in the experiments are larger than Lc, while for the elon- 
gated grains the opposite is true. Since our model has the 
implicit assumption of zero inertia it is inevitable that we 
will only be able to investigate the regime i <C Lc, where 
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SOC is observed. 
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FIG. 5. a) Log-log plot of the probability density of the 
sizes m of the discharge events for several system sizes. It is 
visually apparent that the distribution is bounded, b) Data 
collapse of the curves shown in (a) according to Eq. (y) with 
the exponents k — jj, ~ 1.20. As a visual aid, we display a 
line corresponding to an exponential dependence. 

In summary, we present a new physically motivated 
model for piles of granular material. We find that the 
model self-organizes into a critical state with distribu- 
tions for most quantities described by power laws. We 
measure the exponents characterizing these distributions, 
discuss scaling relations, and find that our model belongs 
to a new universality class. 

We acknowledge discussions with J. Krug, M. Mar- 
kosova, K. Sneppen, H. E. Stanley, and S. Zapperi, 
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